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We present Direct Numerical Simulations of decaying Magnetohydrodynamic (MHD) 
turbulence at low magnetic Reynolds number. The domain considered is bounded by 
periodic boundary conditions in the two directions perpendicular to the magnetic field 
and by two plane Hartmann walls in the third direction. High magnetic fields (Hartmann 
number of up to 896) are considered thanks to a numerical method based on a spectral 
code using the eigenvectors of the dissipation operator. It is found that the decay pro¬ 
ceeds through two phases: first, energy and integral lengthscales vary rapidly during a 
two-dimensionalisation phase extending over about one Hartmann friction time. Dur¬ 
ing this phase, the evolution of the former appears significantly more impeded by the 
presence of walls than that of the latter. Once the large scales are close to quasi-two di¬ 
mensional, the decay results from the competition of a two-dimensional dynamics driven 
by dissipation in the Hartmann boundary layers and the three-dimensional dynamics of 
smaller scales. In the later stages of the decay, three-dimensionality subsists under the 
form of barrel-shaped structures. A purely quasi-two dimensional decay dominated by 
friction in the Hartmann layers is not reached, because of residual dissipation in the bulk. 
However, this dissipation is not generated by the three-dimensionality that subsists, but 
by residual viscous friction due to horizontal velocity gradients. Also, the energy in the 
velocity component aligned with the magnetic field is found to be strongly suppressed, 
as is transport i n this direct ion. This results reproduces the experimental findings of 
iKolesnikov fe Tsinobeil ( 19741) . 


Key words: Low Brn Magnetohydrodynamics, freely decaying turbulence, turbulence di¬ 
mensionality, vortex dynamics, two-dimensional turbulence, quasi-two dimensional flows. 


1. Introduction 

This work concerns the decay of MagnetoHydroDynamic turbulence in electrically 
conducting fluids subjected to an externally imposed magnetic field. We are particu¬ 
larly interested in the influence of solid, electrically insulating walls on this process. This 
generic problem is relevant to a number of practical engineering problems in the metal¬ 
lurgy and nuclear industry, but also bears relevance to some aspects of the dynamics of 
liquid planetary cores and the associated dynamo problem. 

When the magnetic field Be z is imposed (in the sense of the low magnetic Reynolds 
number approximation (Robertd ( 1967I )). turbulence evolves as a result of the competition 
between inertia and the diffusion of momentum along the direction of the magnetic field. 
A structure of size l± becomes elongated over a length l z by this diffusion over a timescale 
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of tj(I z /1±) 2 ( Sommeria fc MorearJ (.1982;)), whilst loosing energy through Joule dissi¬ 
pation (tj = p/(aB 2 ) is the Joule dissipation time, p and a are the fluid density and 
electric conductivity. 1. fMoffatd (1967) first showed that under this linear phenomenology, 
the turbulent kinetic energy decayed at E ~ t -1 / 2 towards an asymptotic state where 
the flow quantities did not vary along the magnetic field (in this sense,a two-dimensional 
state) but where the kinetic energy of the component along the magnetic field was a 
third of the t otal k inetic energy (for a three-component flow). This phenomenology was 
recovered bv 1 Schumann ( 19761) : this author conducted low-resolution direct numerical 
simulations to confirm that this linear phenomenology applied during less than tj but 
that non-linear effects subsequently led to a decay of the kinetic energy associated to the 
velocity along B. However, he also found that the skewness tended to a finite value in the 
later stages of the decay, which he attributed to the persistence of transport of kinetic 


energy B. The more recent and higher resolution simulations of Burattini et al. ( 2010 ) 


confirmed Schumann’s findings. Both studies analysed cases where interaction parameter 
N = tu/tj spanned a range between 0.1 and 50 (tjj(I) = 1/U(Iq) is the eddy turnover 
time based on the initial size and velocity of the large scales Iq an d U(ln)) . Usin g the in¬ 
variance of the ’’parallel” component of Loitsyansky’s integral /||. [Okamoto et all ( 2010h 
showed that for N >> 1, the i -1 / 2 law for the decay of kinetic energy was recovered and 
that the integral lengthscale in the direction of the magnetic field evolved as l z ~ t 1 / 2 . 
At moderate values of N, a similar approach led the authors to conclude that energy 
decayed as E ~ while the integral lengthscales along and across the magnetic field 

increased respectively as l z ~ t 5 / 7 and l± ~ f 3 / 14 . These theoretical scalings as well as 
the invariance of /|| were verified by means of direct numerical simulations at the highest 
resolution available to date (up to 2048 3 ), for N < 1. 

Aside of inertia, a second major factor is likely to interfere with Moffatt’s linear theory: 
the presence of walls, and in particular Hartmann walls, that are perpendicular to the 
magnetic field. These are indeed a feature of practically any of the situations where low 
Rm MHD turbulence is likely to be found. A strictly two-dimensional state is not possible 
in their presence becau se of the very t hin Hartma nn boundary layers that dev elop along 
them (see for instance iMoreaul ( 1990h ). Instead. ISommeria fe Moreau ( 19821 ) theorised 
that in a channel of width L, a structure of size l± became quasi-two-dimensional after 
t 2d(1±) ~ tj{L/1 ± ) 2 . Past this stage, electric current in the core became of order Ha _1 , 
the ratio of the boundary layer thickness to L: dissipation occurs then almost exclu¬ 
sively in the boundary layers and is equally viscous and magnetic. In contrast, strictly 
two-dimensional states are possible when walls are absent and_the Joule dissipation can 
therefore drop to much lower values. Kolesnikov fe Tsinobeil ( 19741 ) also observed exper¬ 
imentally that in the presence of Hartmann walls, transport along the magnetic field was 
suppressed in the later stages of the decay. This was interpreted as an evidence of suppres¬ 
sion of the velocity component in this direction, in contrast to the prediction of theories 
and simulations where no wall was present. Nevertheless, although a ’’through” velocity 
component is precluded by the walls, Ekman pu mping can still potentially lead to strong 
vertical velocities at moderate N ( Potherat et al\ 20001 ). More recently, it was also found 
that even for N > 1, a small amount of three-dimensionality coul d lead to a com plex 
system of three-d imensional co- and contrarotating recirculations (jPotherat et «/Jl201.') : 
Baker et al. [201 5). 


Until recently, numerically simulating MHD turbulence at high N in the presence of walls 
incurred prohibitive computational costs because of the need to resolve the Hartmann 
boundary layers. Recently, the authors took a different approach to the simulation of 
these flows based on spectral methods using bases of functions whose elements already 
























































Figure 1. Geometry of the Channel flow with transverse magnetic field 


incorporate these layers. These partly alleviate this computation al constra i nt to t he point 
where the computa t ional cost becomes independent of Ha ( Dvmkou fc Potheratl ( 20091) : 


Kornet fe Potheratl ( 20141) ). We propose to take advantage of this new technique to in¬ 


vestigate decaying turbulence in a channel bounded by electrically insulating walls in 
view of answering the following questions. 

(a) Does three-dimensionality subsist in the later stages of the decay (t » t 2 d(I_l))? 

( b ) Which part of the energy subsists in the third velocity component ? 

(c) How do Hartmann walls affect the early phases of the decay (t < T2 d(1 ±)) ? 

We shall first recall the governing equations and the timescales that govern the problem 
(section [2j). Our numerical method and simulation strategy is presented in section [3] 
We then examine the earlier decay phase which is expected to present the strongest 
similarities with earlier works not involving walls (section |4|). The later stages of the 
decay where similarities with two-dimensional turbulence are expected are analysed in 
section[5] The robustness of our results is tested in sectionlHlby changing initial conditions 
and domain size. 


2. Governing equations 

2.1. Problem definition 

At low Magnetic Reynolds number, the full system of the induction equation and the 
Navier-Stokes equations for an incompressible fluid can be approximated to the first order 
(The Magnetic Reynolds number Bm represents the ratio of the induced magnetic field 
to the imposed one). This leads to the following system ( Roberts! 1967 ): 


du 1 1. 

— + (u • V)u = —Vp + uAu + -j x B , 
at p p 

V ■ u = 0, 

V-j = 0, 

j = <r(—V $ + u x B), 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 

(2.4) 


where u denotes fluid velocity, B - externally imposed magnetic field, j - electric current 
density, v - kinematic viscosity, a - electrical conductivity, <f> - electric potential. We con¬ 
sider a channel flow with a homogeneous transverse magnetic field Be z and impermeable 
(u|„, Q zz = 0), electrically insulating (j ■ n|,„ a zz = 0) walls located at z = ±L /2 (see fig. Q]). 
In the x y dire ctions we impose periodic boundary conditions with period L. Following 
Roberts! 1 19671) . the Lorentz force can be expressed as the sum of a gradient of magnetic 
pressure p m and a rotational term: 
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j x B = -Vp m - aB 2 A~ 1 d zz u. 


(2.5) 


Using the above identity and adopting the reference scale L , time L 2 /v and velocity v/L 
the set of equations (I2.UI2.4I) can be expressed in dimensionless form: 

^ + P[(u ■ V)u] = Au - , (2.6) 

where Ha = LB^faJpv is the Hartman number and P denotes orthogonal projection 
onto the subspace of solenoidal fields. 


2.2. Timescales of the decay 

We now consider an initially turbulent, isotropic flow, left to decay from t = 0 under 
the action of the Lorentz force and viscous friction, in the configuration described above. 
Hartmann walls can be expected to exert little influence on the initial phase of the 
decay, during whi ch strong two-dim e nsionalisation o f the fl ow should take place, as in 
the simulations of Schumann ( 19761) : lOkamoto et al. ( 2010|j and others. Unlike, in this 
previous studies, the presence of physical walls may lead to a physically realistic phase 
beyond this initial one, where the flow dynamics may become a two-dimensional to some 
extent. To obtain a first estimate for the timescale for the transition between these two 
phases, let us first consider a single turbulent structure of size l±, which diffuses across 
the channel width L under the action of the Lorentz force in time T2d(1±) = tj(L/1±) 2 . 
If this timescale is shorter than both the inertial timescale tu( 1±) = l±/U(l±) and the 
viscous timescale t v {1±) = then this single structure is quasi-two-dimensional. Two 

conditions for two-dimensionality ensue: 


(l±\ 3 lo U(l ±) 1 

\LJ LU(l 0 )N(l 0 y 



(2.7) 

( 2 . 8 ) 


where lo denotes the size of the large scales at t = 0. Since both these conditions are 
scale-dependent, each of them defines a minimum two-dimensional scale. Consequently, 
in a turbulent flow where N = tjj{1±)/tj » 1 , larger scales are two-dimensional while 
smaller scales may still be three-dimensional. From inertial condition m , the small¬ 
est quasi-two-dimensional structure satisfies 1^° ~ L[(l 0 U(l 2 i D )/(N(lo)LU(lo))] 1 ^ 3 ■ Since 
the two-dimensionalisation time T 2 d( 1 _l) increases with l±, this scale is also the slowest 
structure to become two-dimensional and so = T 2 d{Ij °) provides an estimate for the 
two-dimensionalisation phase: 

~ ri { N{la) ^md' ■ (M) 

t 2 D depends on U{1']°)- which is expected to drop by several orders of magnitudes during 
the two-dimensionalisation phase. Not only does this variation of U(l^°) considerably 
slow down two-dinrensionalisation at the initial scale lj 0 (t = 0) , but scales smaller than 
the initial value of lf D may satisfy (12.71) and may in turn become quasi-two-dimensional 
in a time significantly longer than the initial value of (12.91) . On these grounds, a lower 
estimate for the two-dimensionalisation time is obtained by evaluating (12.91) based on the 
value of U(l ±) at t = 0 (this value is fixed by the choice of initial turbulent spectrum). 

From viscous condition (12.81) . by contrast, the smallest quasi-two-dimensional scale 1]° 
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does not depend on U(l±), and neither does the associated two-dimensionalisation time: 

l 2 j° ~ LHa~ 1/2 (2.10) 

= t 2 d{1 2 i D ) ~ rjHa -1 = t h (2.11) 

Since no scale smaller than (12.111) can become two-dimensional regardless of how much 
turbulent intensity drops, th represents a closer estimate of the timescale for two- 
dimensionalisation of the whole turbulent flow than (12.91) . 


2.3. Two-dimensional decay 

Once all structures have become quasi two-dimensional, the evolution of the flow is 
govern ed by two-d i mensio na l dyn amics with an added friction due to the Hartmann 
layers. [Sommeria fc Moreau (1982) showed that the velocity averaged across the channel 
u = /_ 1 u±dz satisfied a shallow water equation of the form: 


<9u 1_ . _ 2u 

— + (u- V ± )u = -~V±p + nA±u - 

at p t h 

V_l • u = 0, 


( 2 . 12 ) 

(2.13) 


where operators with subscript _L operate in the x-y plane only, th appears as the 
typical time for the dissipation due to the Hartmann boundary layers, and therefore a 
characteristic time of the two-dimensional dynamics. From (12.121) . the evolution of the 
total kinetic energy E ~ E 2 D = ||u ||2 associated to the mean flow u when the flow follows 
a two-dimensional dynamics reduces to 


1 dE2D 

2 dt 


= -2 


E2D 

th 


- HlVu" 2 


2D) 


(2.14) 


where || • I 2 d represents the two-dimensional C 2 norm. Introducing lengthscale If = 
/l|Vu|| 2 d ) 1 / 2 which characterises velocity gradients in the (x,y) plane, it comes 


that 


I2D/ 


dE2D 

dt 


4 

th 


1 + 


1 


2 Ha V l 


E2D- 


(2.15) 


It follows from the respective definitions of the total kinetic energy E and E 2 D , that 

E = E 2D (1 + d(max{iJa _1 , a 2 })) , (2.16) 

where a = ||u — u||/||u|| represents the degree of three-dimensionality in the flow. For a 
quasi-two-dimensional flow, E = £?2d(1 + O^Ha -1 )). We shall see from the analysis of 
flow profiles in section [5] (figure [6]) that in the later stages of the decay, a < 0.1 so E 2 d 
can be expected to provide a good approximation for E in the two-dimensional phase of 
the decay. 

The two-dimensional dynami cs o f the flow is expected to favour the formation of large 
scales and indeed ISchumann ( 19761 ) showed that energy transfer towards them occurred 
during the decay. Areas of strong shear may however persist between them. Further¬ 
more, the typical leng thscale of t he vis cous core of quasi-two-dimensional MHD is known 
to scale as LHa ~ 1 / 2 ( Sommerial 1988j). For such fine quasi-two dimensional structures, 
Hartmann friction and horizontal viscous friction would be of the same order. For large 
structures on the other hand, lf/L should be of the order of unity and so in the limit 
Ha —> 00 , the decay should be strongly dominated by Hartmann friction. The total ki¬ 
netic energy should then decay as E ~ exp(—4 t/rn)- Any discrepancy to exponential 
decay of this form is therefore the signature either of thin quasi-two-dimensional struc¬ 
tures or of a residual three-dimensionality. We shall attempt to measure this discrepancy 
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in our numerical simulation to identify the mechanisms of the long-term decay. It should, 
however be noted that for structures such that l±/L ~ ffa -1 / 2 , three-dimensionality 
subsists anyway because at this scale, friction between horizontal planes balances the 
diffusion of momentum along the magnetic field. 


3. Numerical approach 


3.1. Numerical method 


The problem set out in section 12.11 is solved numerically, using a new type of spectral 
method designed to alleviate the computational cost associated with strong anisotropy 
and thin Hartmann boundary layers. Thanks to it, increasing the magnetic held incurs 
essentially no direct computational cost per time step. The mathem atica l found ations of 
this m e thod and numerical implem entations are described in detail in lDvmkou fc Potherat 
( 2009h : lKornet, fc Potheratl ( 2014h . where it is also tested for the exact channel geometry 
studied here. For the sake of completeness, we shall nevertheless outline the principle of 
this new method. Using the spectral approach we seek the solution of eq. m as the 
decomposition on elements of basis U;: 


u = y>(f)m(x). 


(3.1) 


As the spatial dependence is carried solely by m, when representation m is injected 
into eq. mm . the latter reduce s to the set of o rdinary differential equations in time on the 
unknown coefficients Ci(t ) (see lCanuto et al\ ( 20061 ) for a detailed description of spectral 
methods). 

For the basis u; we choose the set of eigenvalues of the operator C formally expressed 
as the right hand side of eq. ( 12 . 61 ) . with the electrical and kinematic boundary condi¬ 
tions of the problem. These functions are a natural choice as elements of a functional 
basis, because the features of flows at high Ha are strongly determined by the prop¬ 
erties of this operator. For example they include specific features of the flow such as 
laminar and turbulent Har t mann boundary layers tha t develop along the channel walls 
( Dvmkou fe Potheratll200ll Potherat fe Dvmkoul 2010l) . Moreover, these modes all have 
negative eigenvalues, and it can be shown that to resolve the flow completely, it is only 
necessary to take into account all modes with eigenvalue A of mo dulu s be lo w a maximum 
|A max |, such that their total number scales as Re 2 /Ha ( Potherat fe Alboussierel 120061) . 
where Re is the Reynolds number based on the large scales. Since the operator C repre¬ 
sents the sum of viscous and ohmic dissipation, the set of modes defined in this way is 
in fact the set of least dissipative modes. 

The main difficulty in solving equation (12.61) using the least dissipative modes lies in 
calculating the spectral representation of non linear terms G(u(xi,yi, Zi)). We use a 
pseudospectral approach and calculate these terms in real space. Therefore we need a 
method to reconstruct the spectral coefficients g n of physical vector fields known at a 
discrete set of points in space x^. To this aim we first use the fact that the eigenmodes 
of C can be factorised as the product of two scalar functions of x and y respectively, and 
a vector function of z. Moreover, the functions of x and y consist of Fourier modes, so 
the set of eigenmodes can be enumerated by a tuple of three numbers (n x ,n y ,n z ) and 
for every mode we can define the vector function En x ,n y ,n z (z) such that each mode takes 
the form 


E. 




(z) exp ( ik nx x + ik n y). 


(3.2) 


























7 


The decay of wall-bounded MHD turbulence at low Rm 

Therefore we first calculate the two-dimensional Fast Fourier transform in the x - y 
directions. This brings the transformed non linear terms under the form: 

G(u(xi, yt, Zi )) = ^2 A %,» s ( z i) ex P (i2nn x Xi + i2nn y yi ), (3.3) 

Tlx •t'ri'y 

where A. nxtUy is the complex amplitude of Fourier mode ( 2Trn x ,2nn y ). Then, for every 
value of ( n x ,n y ) we find the set of spectral coefficients {g nx ,n y ,n z } by solving a set of 
equations 

^ ^ 9n x ,n y ,n z ^n x n y n z (Zj) = A n x ,n y {Zi) . (3-4) 


As the coefficients in this set of equations are constant during a single numerical run, it 
is worth performing the LU decomposition of the corresponding matrix at the beginning 
of calculations and later use it to efficiently find the spectral decompositions. Finally 
the projection onto the subspace of solenoidal vector fields is done by neglecting the 
coefficients corresponding to modes with non zero divergence. Using the Fast Fourier 
transform in x— y planes imposes the distribution of discretisation points in these planes: 
they have to form a regular rectangular grid. We denote its dimensions as N x x N y . In 
our simulations we also use a uniform grid in z direction of dimension N z . For the set of 
equations (13.41) to have a unique solution, the number of modes used during the spectral 
decomposition has to be equal to N z , and the total number of independent modes used 
in the calculations is N = N x N y N z . The technique described above has the advantage 
that the obtained spectral decomposition reproduces exactly the physical field on the 
given set of discretisation points. Therefore momentum and energy are conserved by 
this procedure. However the spectral coefficients g n obtained in this way are different 
from the exact ones 'g ni that would be obtained by decomposition of the same vector 
field over the space of infinite dimension spanned by all eigenvectors of C. \g n — g n \ 
is the so-called aliasing error. To correct this e rror w e adapt the 3/2 technique known 
from standard spectral methods (jCanuto et al. 20061 1. Namely we perform the discrete 
transformation with additional number of modes N larger than the one strictly required 
by the system’s dynamics, Njj (The latter is of the o rder o f the attractor dim ensio n 
of the dynamical system underlying the given problem ( Potherat fc Alboussierell2006l) ). 
After every evaluation of the spectral decomposition, the coefficients corresponding to 
these additional modes are set to 0. 

The spectral method de scribed ab ove was i mplemented by modifying the spectral code 
TARANG developed by Verma et al. ( 20131 ). 


3.2. Simulation strategy 

The bulk of our numerical simulations was based on a domain made of a cube of dimen¬ 
sion L divided uniformly into N x , N y and N z cells respectively in x, y and z directions. 
In order to limit the dealiasing errors we always resolve each of the Hartmann layers 
with at least three computational cells in the 0 direction. Our strategy to study the de¬ 
cay of MHD turbulence relies on four different types of simulations, all gathered in table|T| 
The first type is insp ired from the DNS of decaying MHD turbulence in a three-dimensional 
periodic domain bv [Okamoto et all ( 2010h : the initial conditions consist of a isotropic, 
random Gaussian velocity field with u(k ) ~ exp [(— k/k p ) 2 ] where k p = 47 t/L. This corre¬ 
sponds to an energy spectrum E ~ A: 4 exp [— 2{k/k p ) 2 }. For this choice of initial velocity 
field, the initial integral scale of turbulent motion is given by Iq = \[2/n/k p . The velocity 
spectrum was normalised in such a way that cell sizes in x and y directions correspond to 
Ik/ 1.4 where Ik = lRe~ 3 / 4 is the Kolmogorov length scale and the Reynolds number in 
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its definition Re = u'l/v is based on l and velocity v! = u{k = k p ). This strategy allows 
us to calculate the most intense flow possible whilst minimising mesh-induced numerical 
errors at a given mesh size, since the mesh is always uniform. 

To characterise the influence of the walls, we performed additional simulations starting 
from exactly the same initial conditions as before, but with periodic boundary conditions 
imposed in all three directions. For this set of calculations we used a traditional spectral 
code, Turbo, which uses F ourier modes as the functional basis and was tested and op¬ 
timised for Low -Rm MHD ( Knaepen fc Mom 2004 ; Vorobev et al. 2005 ). 

To evaluate a possible influence of the size of the domain in the x and y dimensions, we 
also performed several simulations with a domain of dimension 2 L in these directions. 
These were computationally expensive and therefore run over a shorter period of time 
than the simulation over a cube of size L. 

Finally, The simple order of magnitude analysis from section 12.21 shows that under a 
strong magnetic field, the decay of MHD turbulence can occurs over a time of the order 
of a few times tj, while a truly quasi-two-dimensional behaviour would not be expected 
before smaller scales are diffused over the height of the channel, i.e. at much later times, 
of the order of tjj = Harj. Thus, if the initial integral scale Iq is chosen much smaller 
than the size of the box, then turbulence will have lost practically all of its kinetic energy 
by the time the two-dimensional dynamics potentially becomes dominant. This would 
make it difficult to study the later phase of the decay. On the other hand, Iq/L needs to 
be sufficiently smaller than unity for a significant three-dimensional phase of the decay to 
exist and so as to generate a sensible transition between three- and two-dimensional dy¬ 
namics. To reconcile these antagonistic constraints, we chose Iq/L = l/^/[2n) ~ 0.4 and 
also run additional simulations initialised with the velocity field obtained at t = 0.5 tj,tj 
in the previous simulations, but in which the velocities are renormalised, so that the total 
energy is restored to its level at t = 0. The flow obtained this way is much closer to a 
quasi-two-dimensional one than the simulations initialised with a random field. Compar¬ 
ing the evolutions of these two different types of initial conditions shall thus give us a 
good measure of the robustness of the features observed in the later stages of the decay 
to a change of initial conditions and to the intensity of the initial flow. 


4. Three dimensional phase 

From the evolution of the total energy (figure [3]), the flow progresses in every run 
through three consecutive phases: during the first one, the flow adjusts from the initial 
conditions. This phase is very short, (shorter than 0.05 tj in all cases). In this section, 
we shall characterise the phase that immediately follows, which features strongly three 
dimensional turbulence, a fast energy decay and lasts several Joule times. For the time 
being, the analysis shall be restricted to cubic domains of size L with Hartmann walls, or 
with periodic boundary conditions when specified. This phase can be identified through 
the strong three-dimensionality visible in the spatial RMS for all (x, y ) of the profile 
along x of magnitude of uj^ on fig. [2] Significant variations along the z direction exist 
until approximately t ~ 1.5 — 2t 2 d{1) — 5 — 10tj for all values of Ha This reflects the 
prominent of the contribution of the large scales to the RMS velocity fluctuations, and 
confirms that large scales indeed become quasi-two-dimensional in this typical time. 


4.1. Total kinetic energy 

Fig. E] presents the evolution of the total kin etic ener g y in t his phase for different values 
of Ha. To compare this evolution to Okamoto et al. (2 0ld )’s laws for the decay of un- 
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Figure 2. Left column: evolution of the normalised, spatial RMS of all vertical profiles of 
< > over all (a;, y) in the domain. Right column: evolution of u z along a vertical line in the 

middle of the domain. 
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Ha 

N x x Ny x N z 

Boundary conditions in z 

N(t = 0) 

o 

II 

f 

Energy Boost 

L X y / Lz 

112 

340 x 340 x 340 

non slip, insulating 

8.55 

336 


N/A 

1 

224 

340 x 340 x 680 

non slip, insulating 

34.2 

336 


N/A 

l 

224 

340 x 340 x 680 

non slip, insulating 

39 

261 

t 

= 0.5tj 

1 

224 

340 x 340 x 680 

non slip, insulating 

43 

229 

t 

= 1 TJ 

1 

448 

340 x 340 x 1340 

non slip, insulating 

137 

336 


N/A 

1 

448 

340 x 340 x 1340 

non slip, insulating 

161 

246 

t 

= 0.5tj 

1 

448 

340 x 340 x 1340 

non slip, insulating 

170 

225 

t 

= 1 TJ 

1 

448 

680 x 680 x 1340 

non slip, insulating 

97 

238 


N/A 

2 

896 

340 x 340 x 1340 

non slip, insulating 

548 

336 


N/A 

1 

896 

340 x 340 x 1340 

non slip, insulating 

661 

230 

t 

= 0.5 Tj 

1 

896 

340 x 340 x 1340 

non slip, insulating 

694 

213 

t 

= Itj 

1 

896 

680 x 680 x 1340 

non slip, insulating 

387 

238 


N/A 

2 

112 

340 x 340 x 340 

periodic 

8.55 

336 


N/A 

1 

224 

340 x 340 x 340 

periodic 

34.2 

336 


N/A 

1 

448 

340 x 340 x 340 

periodic 

137 

336 


N/A 

1 

896 

340 x 340 x 340 

periodic 

548 

336 


N/A 

1 


Table 1. Summary of parameters of calculated 3D cases. 




bounded, three-dimensional and initially isotropic MHD turbulence, we have fitted the 
evolution of energy to laws of the form a(l + bt) c in ranges from 0.05t j to up to 2 tj 
(a, b and c are real constants). All values of c are presented in table [2] lOkamoto et al. 

1 2010h showed that in the limit N —t oo, exponent c should be equal to 1/2. In our cases, 
we obtain the best fits for c = 0.95, 0.75, 0.69 and 0.50 for Ha = 112, 224, 448 and 896 
respectively, over [0.05rj, tj ]. This fit is also relatively robust to a variation of the fitting 
interval, as exponents decrease only slightly when the interval is extended to [0.05tj, 2 tj\. 
From this, we infer that the decay of turbulence between walls is in this line with these 
authors’ prediction over a duration of about tj. Values of c are however slightly higher 
for cases with periodic boundary conditions over this interval, which suggests that the 
influence of the walls is present but moderate for t < tj. This phase corresponds to 
roughly 25-50% of the time interval where we identified strong three-dimensionality in 
the profiles of RMS velocity fluctuations (figure [2]). From times t > tj, by contrast, the 
kinet ic energy tends to decay at a slightly slower rate than predicted bv lOkamoto et al .1 
1 2010h in all cases. This is an indication that some of the turbulent structures interact 
with the walls, as they become stretched vertically under the effect of diffusion by the 
Lorentz force. The energy of such structures is dissipated partly by the action of eddy 
currents recirculating in the Hartmann layers. This Hartmann friction mechanism is typ¬ 
ically Ha times slower than Joule dissipation, which would be the unique electromagnetic 
dissipation mechanism if these structures were not in contact with the wall. This explains 
that the decay of energy slows down for t>Tj. 

All cases where the energy was boosted during the three-dimensional phase exhibit a 
similar behaviour to cases where the energy was left to decay from the start. Neverthe¬ 
less, for all of them, the value of the c coefficient fitted over intervals of tj or more is 
lower than for their counterpart without energy boost. Furthermore, the later the en¬ 
ergy is boosted, the lower the value of c. This slower decay reflects the influence of the 
anisotropy of the boundary conditions: at the time of the energy boost, the flow recov¬ 
ers the same energy as the initial one but conserves the anisotropy that has developed 
during the initial decay, before the energy was boosted. Consequently, vortices are more 
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Figure 3. Evolution of total kinetic energy, normalised by its value at t = 0 (solid line). The 
dashed line represents the fitted law of the form a(l + bt) c over interval [0.2, 2 tj]. 


elongated, interaction with the walls is more significant and the slower friction in the 
Hartmann layers represents a more important fraction of the dissipation. 

4.2. Dissipation 

Initially, the kinetic energy is mainly dissipated ohmically and the initial ratio of viscous 
to Joule dissipation scales as ~ 1/Ha. As the flow becomes more two-dimensional, both 
viscous and Joule dissipations diminish in the bulk of the flow (between the Hartmann 
layers). Conversely, dissipation in the Hartmann boundary layers increases as more and 
more structures interact with the walls. Therefore, the main contribution from the total 
dissipation ends up coming from the Hartmann layers (see fig. []}. For Ha = 112, 224, 
448 and 896 dissipation in the layers becomes larger than in the bulk at t = 25 tj, 37.6 tj, 
45.1tj and 54tj respectively. However in the Hartmann layers the Joule dissipation is 
nearly the same as viscous dissipation. Therefore the global ratio of viscous to Joule 
dissipation increases with time. The value of this ratio becomes larger than unity for 
t = 6.3 tj, 19 tj, 50.1 tj and 88tj (or equivalently at t = 0.0567#, 0.084T#, O.llr# and 
0.27#) for Ha = 112, 224, 448 and 896 respectively. In all cases, the ratio e„/ej tends to 
an asymptotic value of ~ 1.3 after a time of the order of 7#. From this perspective, the 
dissipation behaves as in a three-dimensional flow during a period of time that is longer 
than the timescale of three-dimensional Joule dissipation tj. but shorter than that of two- 
dimensional effects r#. This intermediate time scale is an indication that although large 
scales become two-dimensional over a time of the order of T 2 d( 1 ), smaller scales remain 
three-dimensional during a significantly longer period of time. After t > 1.5 — 272#, the 
contribution of the large scales to the total dissipation comes mostly from the Hartmann 
layer, where it is weak, whereas the contribution of the small scales consists of stronger 
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E fitted from 0.05tj to I l z fitted from 0.05tj to 


Ha Remarks 

2 TJ 

TJ 

0.5tj 

0.2tj 

2 TJ 

TJ 

0.5tj 

0.2tj 

112 

-0.91 

-0.95 

-0.97 

-0.86 

0.18 

0.21 

0.24 

0.26 

224 

-0.65 

-0.75 

-0.82 

-0.71 

0.26 

0.33 

0.37 

0.42 

448 

-0.59 

-0.69 

-0.42 

-0.83 

0.27 

0.34 

0.38 

0.44 

896 

-0.51 

-0.51 

-0.41 

-0.34 

0.26 

0.34 

0.37 

0.44 

224 E boosted at 0.5 tj 

-0.54 

-0.58 

-0.63 

-0.68 

0.20 

0.25 

0.28 

0.29 

448 E boosted at 0.5rj 

-0.46 

-0.53 

-0.73 

unstable 

0.20 

0.26 

0.29 

0.40 

896 E boosted at 0.5 tj 

-0.45 

-0.56 

unstable 

unstable 

0.21 

0.29 

unstable 

unstable 

224 E boosted at 1 tj 

-0.50 

-0.51 

-0.52 

-0.56 

0.15 

0.20 

0.22 

0.22 

448 E boosted at 1 tj 

-0.43 

-0.48 

-0.70 

unstable 

0.18 

0.23 

0.25 

0.39 

896 E boosted at 1 tj 

-0.43 

-0.52 

unstable 

unstable 

0.19 

0.24 

0.36 

unstable 

112 periodic BC 

-0.76 

-0.86 

-0.96 

-1.05 

0.33 

0.37 

0.40 

0.42 

224 periodic BC 

-0.56 

-0.69 

-0.79 

-0.74 

0.40 

0.50 

0.55 

0.59 

448 periodic BC 

-0.52 

-0.65 

-0.82 

-1.02 

0.42 

0.50 

0.56 

0.64 

896 periodic BC 

-0.47 

-0.58 

-0.66 

-0.89 

0.43 

0.51 

0.57 

0.65 

448 Domain 2L x 2L x h 

-0.91 

-1.04 

-1.14 

-1.22 

0.38 

0.49 

0.56 

0.66 

896 Domain 2L x 2L x h 

-0.87 

-1.03 

-1.29 

-1.68 

0.37 

0.46 

0.45 

0.38 


Table 2. Fitted values for constant c in lOkamoto et all (|2010T )’s laws of the form a(bt+ l) c for 
the decay of total kinetic energy and the growth of integral scale l z . Cases where the formal 
error of the least squares method was greater than 20% are marked as ’’unstable”. 


Joule dissipation in the bulk. Joule dissipation therefore remains higher than viscous 
dissipation until this contribution has significantly reduced, which occurs on a timescale 
of at most th (see section & 


4.3. Integral lengthscales 

Figure [5] (top) shows the initial evolution of the integral lengthscale in the 2 direction l z , 
and in the direction orthogonal to the magnetic field lj_. These are respectively defined 
as: 


Lh z (x, y, z)u z {x, y,z + z') d Vdz' 
f u z( x , y, z) dV 

1 (f_h X (x, y, z) ■ u x (x + x', y, z) d Vdx' 

2 \ fu%(x,y,z) dV 

I f u y( x , V , z ) • u y (x, y + y',z) d Vdy' \ 
ful(x,y,z)dV )’ 


(4.1) 


(4.2) 


The initial growth of l z can be again fitted with the formula_a(l + bt) c . Except for 
t < 0.2 tj, the fitted value of c is significantly smaller than lOkamoto et al. ( 2010f) ’s 
theoretical value of 0.5. It is also strongly dependent on the fitting interval for Ha = 224, 
448 and 896. For Ha = 112 the integral lengthscale l z grows even more slowly with a 
fitted value down to c ~ 0.18 over [0.05, 2tj). This behaviour indicates a very strong 
influence of the walls on the growth of l z from the outset of the decay. It is somewhat 
remarkable that despite this early influence of the walls on l z . the energy decay shows 
little influence of the walls during as long as tj - 2 tj. During this initial stage l z also 
grows faster for larger values of Ha , with indication that at Ha = 224, it is already close 
to its asymptotic behaviour (in the sense of large Ha). For t > 1.5rj, the growth of l z 
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Figure 4. Left column: evolution of the kinetic energy in cases with walls (solid lines), periodic 
boundary conditions (dashed line). Right column: ratio of viscous to Joule dissipation in cases 
with walls (solid line) and periodic boundary conditions (dashed line). The dotted line represent 
the fraction of energy which is dissipated in Hartmann layers in cases with wall. On all graphs, 
the shorter curves represent the case with walls where energy was boosted at t = 0.5 tj (Short 
dashed line: Joule to viscous dissipation, short dash-dotted line: fraction of dissipation in the 
Hartmann layers) . 
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t/rj 


(left) and l± (right) in the presence of 


t/rj 

Figure 5. Evolution of integral lengthscales l, 

Hartmann walls. 


slows significantly. As for the energy, this is due to the increasingly wide range of scales 
at which structures reach the walls during the two-dimensionalisation process, and whose 
growth in thus impeded in the z direction. 

With periodic boundary conditions, the fitted value of exponent c is higher than with 
Hartmann walls. It is close to the theoretical value of 0.5 at high Ha for t < tj, and 
decreases thereafter. Unlike in cases with walls, this behavi our is qu ite inse nsitive to the 


fitting interval for t < tj, which confirms the validity to Okamoto et al. ( 2010l) ’s law 


at high Ha over at most one Joule time. The validity of this law also indicates that 
the parameter l±(t = 0 )/L was chosen sufficiently small to observe the main features of 
three-dimensional unbounded turbulence in a periodic domain. However, the fact that 
this exponent is higher with periodic boundary conditions than with walls suggests that 
eddy currents circulating between the Hartmann layer and the bulk (when walls are 
present) strongly increase the influence of the boundaries compared to the periodic case 
where this effect is absent. 

The integral lengthscale l± grow s very slowly during the decay. This is consistent with 
the prediction of lOkamoto et al\ ( 2010t ) of a growth as (t/rj) 1 / 7 . However such a small 
exponent is difficult to quantify on a timescale of the order of tj, where it is expected 
to be valid. Furthermore, l± evolves slowly all the way through our calculations, with 
no clear evidence of different behaviour when the flow is close to two-dimensional than 
when it is three-dimensional. 


5. Quasi-two dimensional phase 

We shall now describe the later stage of the flow evolution where it approaches a 
quasi-two dimensional behaviour, and characterise this asymptotic regime. 

5.1. Velocity profiles 

When this stage is reached, the spatial RMS over all (x,y) of the profile along z of 
uj^ have already been considerably smoothed out during the three-dimensional phase of 
the decay. In every case, we were able to identify a time tq 2 Dj from which the profile 
starts to flatten monotonically without qualitatively changing shape. This type of decay 
would be expected from a flow governed by mostly two-dimensional dynamics. tq 2 d was 
defined as the time at which the maximum value in the velocity profile starts decreasing 
monotonically. It was found at 19.3rj(= 0.34 th), 37.6tj(= 0.336t#), 65.2tj(= 0.29 th) 
and 98tj(= 0.22 th) for Ha = 112, 224, 448 and 896 respectively. The fact that it obeys 
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t/rj 


Figure 6. Evolution of the maximum of < > xy /< > xy (as represented on figure [2}. 


a timescale of about 0.3 tjj that is commensurate with the two-dimensional timescale th 
indicates that the large scales only acquire a two-dimensional dynamics once a significant 
part of the spectrum is close to being two-dimensional (since two-dimensionalisation 
of the whole spectrum is expected to occur over a period of approximately th)- For 
t > tq 2 d, the evolution of the shape of the profiles is practically unaffected by smaller 
scales that still retain a three-dimensional behaviour at this stage of the decay. 
Remarkably, in none of the cases, did we find that the profile was quasi-two-dimensional 
after t > th (although such very long times could not be reached for Ha = 448). Instead, 
all profiles seem to have reached a barrel-like shape after t ~ 50tj, which evolves only 
very slowly after this time. Also the sha pe o f t he pro file is f latter for the larger values 
of Ha. This shape was first theorised bv lPotherat et al\ (2 000 1 and numerically observed 
bv lMiick et al. ( 2000 ). It stems from eddy currents recirculating between the Hartmann 
layers and the core. These currents are driven along the axis of columnar vortices. Their 
leak into the core drives differential rotation betwee n horizontal planes of the vortex 
leading to the barrel-shaped profile (see also jPotheratl ( 2012l l). Since the effect is driven 
by currents recirculating between the bulk and the Hartmann layer, which scales as Ha — , 
it is less marked at high Ha. 

The intensity of the barrel effect can be measured through the relative value of the 
maximum in the profile. The evolution of this quantity is shown on figure [G] First, this 
graph confirms quantitatively that the barrel effect is less pronounced at higher Ha. 
Second, the graph also confirms the invariance of the barrel shape of the large scales 
beyond t ~ 50r,/. This appears to be verified over at least 100 tj, although the total 
energy drops by a factor of up to a couple of orders of magnitude during this interval, 
depending on the value of Ha. This remarkable feature shows that even at high Ha , a 
form of three-dimensionality subsists a large times, even in the large scales. 

Finally, we verified that cases where the energy was boosted at t = 0.5 tj and t = 
tj exhibit the same behaviour, which indicates robustness to initial conditions of this 


scenario. 


5.2. Energy 

In all calculated cases (with and without energy boost), the evolution of the total kinetic 
energy in the late phase can be described as an exponential decay with slowly changing 
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t/rj t/rj 

Figure 7. Evolution of E z /E for cases with walls (left) and with periodic boundaries (right). 


timescale. However even at the latest time in our calculations this timescale does not seem 
to converge towards a clear asymptotic value. This indicates that even for t > tjj, when 
all structures could theoretically be expected to have been diffused across the channel 
height, the decay is still not entirely dominated by dissipation in the Hartmann layers. 
The origin of the extra dissipation shall be determined by examining different types of 
dissipation in section PTTTT1 


The energy associated to the z component of the velocity is represented on figure [7] For 
all cases with Hartmann walls, E z /E tends rapidl y to zero, which is consistent with the 
experimental findings of lKolesnikov fe Tsinoben (119741) for turbulence in a duct. By con¬ 


trast , experiments where turbulence was kept far from Hartmann walls (lAlemanv et al 
119791 1 and simulations with periodic boundary conditions rather than walls (Schumann 


Burattini et al. (12010 )) show a transfer of energy to the velocity component along 


B, which results in a maximum in E z /E , followed by a much sl ower decay t han in the 
case with walls. This supports the explanation put forward bv iBurattini et ~gj\ ( 2010l ) 
who hypothesised that the difference between these two scenarii was due to the presence 
of the Hartmann walls. 

The question of how much transport along B remains asymptotically can be further 
analysed through the evolution of the skewness coefficient 




where 


T = V 2fc 2 q(k) 


k 


*(k), 


k 


2fc 2 u(k) • u*(k), 


(5.1) 


(5.2) 

(5.3) 


which Schumann ( 19761 ). found to remain constant at large times. This author interpreted 
this behaviour as an evidence that transport of E z became important at large times. Our 
computed evolution of S is sho wn on figure [5] In the domain with periodic boundary 
conditions, we recover Schumann ( 1976t l’s findings at high Ha that S seems to converge 
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t/rj t/rj 

Figure 8. Evolution of Skewness S for cases with walls (left) and with periodic boundaries 

(right). 


to a constant asymptotic value. However, as in this authors’ work, we were not able 
to compute S for t > r# at the highest value of Ha, so the two-dimensionalisation 
process may not be entirely complete at the end of this particular calculation. This 
leaves room for the possibility that S may in fact evolve on a timescale significantly 
longer than the duration of our calculations in this case. Simulations with Hartmann 
walls, by contrast, exhibit a fast decay of S and it is readily visible that S —> 0 in all 
calculated cases. Together with the fast decay of E z /E, this behaviour confirms that the 
presence of walls results in a very strong suppression of transport along B, in contrast 
to the phenomenology of flows with periodic boundary conditions. 

Note that in the two-dimensional phase, both E (see fig. HI) and E z /E appear to decay 
faster at lower values of Ha in the case with walls. It should however be pointed out that 
in this phase, turbulence decays over a timescale of r# = Harj, which appears slower at 
higher Ha in units of time normalised by tj when it is in fact faster in dimensional time 
units. 

Finally, in simulations where the total energy in the flow was artificially boosted at 
t = 0.5tj (and, we verified, t = tj) both E and E z /E exhibit a similar behaviour to 
the case without energy boost, which confirms the robustness of this phenomenology to 
initial conditions. 


5.3. Dissipation 

Asymptotically, the ratio of viscous to ohmic dissipations e^/ej tends to a value of 
approximately 1.3 for all values of Ha (see fig. 0. For all cases except Ha = 448 we 
observe a maximum in the temporal evolution of this ratio at times 0.35r#, 0.387"#, 0.5r# 
for Ha = 112, 224 and 448 respectively. The maximum becomes less pronounced with 
increasing Ha, with the ratio of maximal value oi e v /ej to its asymptotic value decreasing 
from 1.28 at Ha = 112 to 1.04 at Ha = 448. We presume that for Ha = 896, this 
maximum is still present for t > 0.5r#, although even less prominent. This behaviour 
is very different from that observed in simulations with periodic boundary conditions, 
where t v jej increases indefinitely. In both cases the sharp initial increase of this ratio 
is due to the two-dimensionalisation of smaller and smaller structures. Indeed, from the 
right hand side of (12.61) . the ratio of dissipations for a structure of wavenumber (k±, k z ) is 
e u (k±, K z )/ej(k±, k z ) = Ha~ 2 ( 1 + ( kj_/n z ) 2 ). While initially k±/n z ~ 1, k±/K z becomes 
very large for nearly two-dimensional structures. As time progresses, this becomes true 
for increasingly large values of k±, and so the ratio of total viscous to Joule dissipations 
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£vltj increases. With periodic boundary conditions, strictly two-dimensional structures 
( k z = 0) can exist and so this ratio is unbounded. In the presence of Hartmann walls, 
on the other hand, quasi-two-dimensional structures generate little dissipation in the 
bulk. Most of their dissipation comes from the Hartmann layers, where viscous and Joule 
dissipation are locally of the same order. This explains that the ratio e^/ej converges to 
a value of the order of unity. 

The presence of the maximum in the evolution of e„/ej a t lower values of Ha may be 
attributed to structures that are too small for magnetic diffusion to stretch them up to 
the walls, i.e. for which at t = 0, (k±) = kj 1 N{k±) 1 / 2 < L. These initially generate 

a large contribution to the ratio e„/ej over a time scale of the order of T 2 d(A_l), which 
is longer than the two-dimensionalisation time of the large scales T 2 d(Iq) but shorter 
than the decay time of quasi-two-dimensional structures th ■ For t ~ th, these structures 
have lost most of their energy and the surviving structures extend across the whole 
channel. As Ha increases, structures that remain three-dimensional during their entire 
life span are confined to a region of the spectrum of higher and higher lengthscale and 
their contribution to the total dissipation progressively vanishes. This explains that the 
maximum is less pronounced at higher values of Ha, and also that it take place at later 
times. 

More details on the asymptotic state can be obtained by inspecting the ratio of the total 
dissipation coming from the Hartmann layers to that coming from the bulk (figure U). 
Asymptotically, the dissipation comes dominantly from Hartmann layers because the flow 
becomes close to quasi-two dimensional at all values of Ha investigated here. However, 
approximately 20% of the dissipation still takes place in the bulk asymptotically. It is 
tempting to check whether the residual three-dimensionality due to the Barrel effect noted 
on the velocity profiles (section 15.11) is responsible for this extra dissipation: since it is 
driven by currents in the Hartmann layers that recirculate in the bulk, Joule dissipation 
must be associated to it. Current conservation implies that the current densities in the 
bulk J c and the Hartmann layers J H must satisfy J H / J c ~ Ha and so the contributions 
to Joule dissipation from the core and the Hartmann layers must be in a ratio of e^/ej ~ 
Ha. This scaling is insufficient to explain the relatively high level of dissipation observed in 
the bulk: this discards the barrel effect as the residual source of dissipation. To find out the 
origin of this extra dissipation, we calculated non-dimensional quantity 2Ha~ l (L/l v A ) 2 , 
which, from (12.151) . represents the actual ratio of quasi-two-dimensional viscous friction 
to Hartmann friction. It turns out that this ratio converges towards a constant value 
around 0.2 regardless of the case considered (see figure flOl) . This implies that extra 
dissipation in the bulk results from viscous friction incurred by gradients of the x and y 
components of the velocity in the x and y directions. This extra dissipation also explains 
that the energy does not decay exponentially even for t > th, as would be expected 
if Hartmann friction was the exclusive dissipation mechanisms. Interestingly, the ratio 
2 Ha -1 (L /1 1 ^) 2 behaves in roughly the same way whether Hartmann walls are present or 
not. Even though Hartmann friction is absent when boundary conditions are periodic 
in all three spatial directions, this ratio still gives a normalised measure of the level of 
two-dimensional viscous dissipation, which appears to be the same with and without 
Hartmann walls. Furthermore, our simulations show that in the Hartmann layers, Joule 
and viscous dissipations are asymptotically the same up to few percent. Therefore viscous 
dissipation due to horizontal velocity gradients also explains that the ratio e^/ej remains 
greater than unity even for t > th, when two-dimensionalisation is expected be complete 
over the whole spectrum of lengthscales. 
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Figure 9. Evolution of interaction parameter N = Bl±/u, 



t/rj 



t/rj 


Figure 10. Ratio of effective dissipation by two-dimensional viscous friction to Hartmann fric¬ 
tion: case with Hartmann wall (left) and case with periodic boundary conditions (right). With 
periodic conditions, this quantity represents the normalised two-dimensional dissipation only, as 
Hartmann friction is absent. 


5.4. Integral lengthscales 

Figure HT1 shows the evolution of l z up to the later stages where t > th- For large values 
of Ha , l z tends asymptotically to a value of approximately 0.8. By contrast, in simula¬ 
tions with periodic boundary conditions l z asymptotically tends to 1. At high Ha the 
Hartmann layers are laminar and not affected by inertia present in the core. The classical 
theory for these layers then implies that the velocity normal to the wall in the layer is 
OlyHa- 1 ) and this imposes a region of very low values of u z near the walls. By contrast, 
residual velocity may exist in the core and so the z-component of the velocity cannot be 
correlated over the entire width of the channel (see profiles of vertical velocity on figure 
0. This effect is absent with periodic boundary conditions where a constant through-flow 
can exist in the z-direction, which allows values of l z close to 1. 

The convergence of l z for cases bounded by walls becomes less smooth for lower values of 
Ha. For Ha = 112, l z shows no signs of convergence to a finite value. It starts decreasing 
around t ~ 45tj and still decreases at f ~ 3 tr- This behaviour is caused by Ekman 
pumping in large quasi-two dimensional vortices. This effect tends to produce vertical 
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t/rj t/rj 

Figure 11. Evolution of integral lengthscales l z (left) and l± (right) for the simulation with 

walls. 



t/rj 


Figure 12. Evolution of integral lengthscale l z for periodic cases. 


profiles of u z that are antisymmetric with respect to the midplane. According to (14.21) . 
reduces the value of l z as its definition is based on u z (z ) (see fig. 0. The i ntensity o f 
Ekma n pumping is driven by inertia but damped by the Lorentz force: IPotherat et all 
1 2000l) showed that it scaled as u z /u± ~ Ha~ 2 N~ l . This explains that this effect is only 
noticeable at the lowest value of Ha. This also suggests that ultimately, since the inter¬ 
action parameter N(t) diverges asymptotically (see figure 0, Ekman pumping should 
progressively disappear and l z may increase again. The evolution of l z on figure [TT] how¬ 
ever implies that this may only take place after a very long time, beyond the reach of our 
calculations. Also, since E z /E is very small in the later stages of the decay, l z relies in 
fact on low values of u z and only reflects a weak component of the flow, from the Energy 
point of view. Despite its relative weakness, Ekman pumping is responsible for the larger 
values of l± at Ha = 224 and more noticeably at Ha = 112 (figure [Till . Ekman pumping 
indeed transports momentum radiall y outward of large structure thus increasing their 
effective size in the horizontal plane ( Sommerialll988l) . 

From this phenomenology, l z appears dominated by the dynamics of secondary flows 
but does not reflect accurately the dimensionality of turbulence. We argue that this calls 
for a more suitable quantity to characterise the growth of vertical scales along the z 
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Figure 13. Evolution of integral lengthscale if for simulations with walls. 


direction. We propose that such a quantity may be defined as: 

7 _l S f u - l(x, y, 2)u_l0t, y,z + z') dVdz' 

fu 2 ± (x,y,z)dV ■ [bA> 

The temporal evolution of if is shown on fig.[l3]for decaying turbulence between walls. In 
contrast to l z it asymptotically converges to 1 for all values of Ha and curves remain very 
close to each other in all stages of the evolution. During the initial stage, if grows over 
a time scale of T2 d( 0> which does reflect the two-dimensionalisation of the large scales. 
if reaches a value of 0.9 at ~ 15 tj for all values of Ha, and its subsequent evolution in 
the quasi-two-dimensional phase is only slow. This indicates if is only weakly affected 
by the two-dimensionalisation of smaller scales, unlike the total kinetic energy and the 
dissipation. 


6. Robustness analysis 

In order to focus on the long term evolution of the decay of turbulence and still 
keep computational costs reasonable, we have conducted our study with a domain of 
limited size and used only one set of isotropic initial conditi ons. Ongoing experiment s 
on decaying turbulence conducted on the FLOWCUBE setup (IPotherat fe Klein (|2014|)) 
show that reliable quantitative laws require ensemble averaging on a large number of 
initial conditions, which cannot be done numerically. Nevertheless, we shall now estimate 
the impact of these choices on the mechanisms found, by examining the result of two 
simulations in a channel four times bigger (dimensions 2L x 2L x L) , with different random 
initial conditions (albeit with the same statistical properties as in the cases where energy 
has not been boosted), and slightly smaller initial Reynolds number (see table [[]). 

Comparison between energy decay and dissipation ratios in small and large domains 
is shown on figure [Td] The decay follows a similar profile in both cases. In the three- 
dimensional phase, fits to Okamoto et al\ ( 2010h ’s decay laws yields exponents that are 
consistently higher in the case of a larger domain. Nevertheless their variations with Ha 
and with the fitting interval are consistent between cases (see table [2J. We also noticed 
that actual values of the exponent are sensitive to the lower bound of the fitting interval, 
but that again, variations with Ha and with the fitting interval are consistent between 
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Figure 14. Evolution of energy (left) and dissipation rates (right) in smaller and larger 
domains (marked ’’big” in the legend), with Hartmann walls. 



Figure 15. Evolution of l z in smaller and larger domains (marked ’’big” in the legend) with 

Hartmann walls 


cases, for the same lower bound. This sensitivity can most likely be attributed to the 
lack of an ensemble ave rage which, as in expe riments, would be required for a precise 
estimate of exponents in Okamoto et al. ( 2010h ’s decay laws. 

Nevertheless, the interval of validity of the these laws, can still be estimated by varying 
the upper bound of the fitting interval, and this yields consistent results between cases 
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and choices of fitting interval. Most importantly, the main properties of the decay outlined 
throughout the paper appear consistent between the two sets of simulations: the duration 
of the three-dimensional phase is of the order of T 2 d(Iq) in both cases, and the asymptotic 
behaviours are identical as decays of energy are parallel to each other. Similarly, Ratios of 
dissipations (dissipation in the Hartmann layer to dissipation in the bulk and viscous to 
Joule dissipation), depict the same phenomenology as the evolution of energy: timescales 
are identical for both domains and evolution curves are parallel. The integral lengthscale 
Ifigure lTbl) too evolves initially in a similar fashion in both cases and converges to the same 
asymptotic value. In conclusion to this short analysis, although i t is diffi cul t to p recisely 
verify the numerical value of the exponents predicted in lOkamoto et al. ( 2010l) ’s laws, 
the scenario for the decay outlined through the analysis cubic domain seems robust to 
changes in numerical parameters. 


7. Conclusion 

Using a new type of spectral methods based on the least dissipative eigenmodes of 
the dissipation operator, we were able to perform direct numerical simulations of freely 
decaying turbulence in a strong magnetic field, between two Hartmann walls. The decay 
exhibits three- and a two-dimensional phases with an overlap: the former is dominated 
by the two-dimensionalisation process, where diffusion by the Lorentz force stretches 
vortices until they reach the Hartmann walls. This process is highly dissipative and leads 
to a rapid variation of energy and of the integral lengthscale along B. Larger scales are 
two-dimensionalised more quickly than smaller ones. Once the large scales of turbulence 
are close to two-dimensional (after T 2 d(Iq )) the flow starts exhibiting a two-dimensional 
dynamics, where dissipation mostly takes place in the Hartmann boundary layers, with a 
slower characteristic timescale th = Harj. However since it can take up to th for small 
scales to adopt a two-dimensional dynamics, there is no clear separation between these 
two phases and both two and three-dimensional dissipation mechanisms co-exist long af¬ 
ter T 2 d(Iq)- We were able to single out several important features of this phenomenology: 
First, the presence of the walls turned out to impede the growth of l z right from the 
earliest stages o f the decay, whereas the decay of energy remained roughly in line with 


Okamoto et al. ( 2010l) ’s law of E ~ t 1 / 2 for unbounded turbulence in the limit of high 


Ha, during around one Joule time. 

Second, energy associated to the velocity component across the channel is very strongly 
suppressed: E z /E tends to 0 much fas t er tha n for unbounded turbulenc e. This result is 
consis tent with Kolesnikov fc Tsinobeil (1974 )’s experiments and confirms Burattini et al\ 
( 2010h ’s hypothesis that the presence of walls is responsible for the suppression of the 
third component. Further evidence of this suppression is visible in the long term be¬ 
haviour of the skewness which tends to 0 in the case with walls. With periodic boundary 
conditions, by contrast, the Skewness apparently tends to a constant value. However, this 
only seems true at high Ha. Since, two-dimensionalisation occurs over a timescale of the 
order of th which is much longer than our calculations at high Ha, and than calculations 
in previous studies, it is unclear whether the Skewness indeed remains constant past the 
two-dimensional phase with periodic boundaries. 

Long into the ’’two-dimensional phase”, we found that even at the highest value of Ha, a 
form of three-dimensionality subsisted, due to currents recirculating between the Hart¬ 
mann layers and the bulk. This effe ct is characteri s ed b y the barrel shape visible on 
the larger structures, as predicted bv IPotherat et al. ( 200011 . Though less pronounced at 
higher values of Ha, our simulations show no evidence of it vanishing at larger times. 
Thirdly, in quasi-two dimensional flows dominated by dissipation in the Hartmann bound- 
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ary layers, the total kinetic energy would be expected to decay exponentially with a 
timescale of th- However, a true exponential decay of this sort was never observed, even 
for t > th- Remarkably this discrepancy to a pure exponential decay did not result from 
the residual three-dimensionality due to the barrel effect, but mostly from viscous friction 
in the horizontal plane. 

Finally, at more moderate values of Ha (Ha = 112), secondary flows in large structures 
significantly affect the decay by increasing the integral lengthscale in the directions along 
the channel (perpendicular to B). Since, however this effect is driven by two-dimensional 
inertia, it is expected to vanish at larger times, but was still present after 2.5 th- When 
present, it is shown to dominate the behaviour of the integral lengthscale in the direction 
of the magnetic field. This prompted us to put forward an alternative definition for this 
integral lengthscale that gives a better measure of the flow dimensionality. 

The authors are grateful to the Leverhulme Trust, who supported this work through 
Research Project Grant ref. F00/732J and to Professor Mahendra K. Verma, who made 
his Fourrier-base spectral code available to them as a basis for the implementation of 
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